# ------------------------------------- #
# Replication code for:
#
# Rathbun, Brian C. and Caleb Pomeroy "See No Evil, Speak No Evil? Morality, Evolutionary Psychology, 
# and the Nature of International Relations," International Organization.
#
# This script reproduces the China survey results from the Appendix.
# ------------------------------------- #


# --- set your working directory --- #
setwd("~/Downloads/seeing_evil_replication/")

# --- Libraries --- #
library(ggplot2)
library(mediation)
library(texreg)
library(plyr)


# --- Beijing Area Study --- #
# columns list characteristics that respondents associate with the people of a given country. 
# the final column "views_of" contains the country about which subjects were being asked.
# characteristics were measured on 7-point scales, with "1" indicating the strongest beliefs about the "good" version of 
# that characteristic, and "7" indicating the strongest beliefs about the "bad" version of that characteristic.
# take the first row, for example: the "views_of" column indicates that the subject was being asked about 
# fellow Chinese individuals. in the "arrogant_modest" column, a value of "1" indicates the subject believes Chinese individuals 
# are "modest." a value of "7" would indicate the subject believes Chinese individuals are "arrogant."
china_bas <- as.data.frame(read_csv("data/china_bas.csv"))

cor_df_temp <- data.frame()
for(i in 1:length(unique(china_bas$views_of))){

  sub_temp <- subset(china_bas, views_of == unique(china_bas$views_of)[i])
  cor_vars <- c("immoral_moral", "arrogant_modest", "insincere_sincere", "barbaric_civilized")
  
  for(k in 1:length(cor_vars)){
    cor_temp <- cor.test(sub_temp$warlike_peaceful, sub_temp[,cor_vars[k]]) 
    
    cor_df_temp <-
      rbind(cor_df_temp,
            data.frame(rho = cor_temp$estimate,
                       ci_high = cor_temp$conf.int[2],
                       ci_low = cor_temp$conf.int[1],
                       variable = cor_vars[k],
                       views_warlike = unique(china_bas$views_of)[i]))
  }
}

cor_df_temp$variable <- factor(cor_df_temp$variable, levels = c("barbaric_civilized", "insincere_sincere", "arrogant_modest", "immoral_moral"))
cor_df_temp$views_warlike <- factor(cor_df_temp$views_warlike, levels = c("china", "us", "japan"))
levels(cor_df_temp$views_warlike) <- c("China\nWarlike", "US\nWarlike", "Japan\nWarlike")
levels(cor_df_temp$variable) <- c("Barbaric", "Insincere", "Arrogant", "Immoral")

# --- Figure A6(A) --- #
ggplot(cor_df_temp, aes(y = rho, x = variable, fill = views_warlike)) +
  theme_bw() +
  facet_grid(cols = vars(views_warlike)) +
  geom_bar(stat = "identity", width = .4, colour="black") +
  geom_errorbar(aes(ymin = ci_low, ymax = ci_high), width = .1, size = .9) +
  scale_fill_manual(values = c("gray40", "rosybrown2", "lightcyan3")) +
  scale_y_continuous(breaks = c(0, .25, .50, .75, 1)) +
  coord_flip() +
  theme(axis.text.y=element_text(size=16),
        axis.text.x=element_text(size=12),
        strip.background = element_rect(fill="gray95"),
        strip.text = element_text(size=16),
        legend.position = "none",
        axis.title.y = element_blank(),
        axis.title.x = element_blank(), 
        plot.title = element_text(hjust = -.1, size = 26),
        panel.spacing = unit(1.2, "lines")) 
#ggsave("bas_correlations.pdf", height = 3, width=7.4)



# --- China USCSPP Survey --- #
# the various "us_[some security issue]" columns contain the importance respondents attach to each
# security issue, with "1" representing "very serious problem" or "somewhat serious problem," 
# and 0 otherwise (indicating the security issue is not important to the respondent).
# "neg morals" column contains the number of immoral traits subjects associate with Americans.
# the remaining columns contain demographic information.
china_moral_threat_df <- as.data.frame(read_csv("data/china_mass.csv"))


dvs <- c("us_arms_sales", "us_mil_pres", "us_suppression", "us_supremacy", "us_tibet", "us_recon")
threat_interaction_plot <- data.frame()
for(i in 1:length(dvs)){
  
  m_temp <-
  glm(eval(as.name(dvs[i])) ~ us_greatest_threat + neg_morals + us_greatest_threat*neg_morals,
      data = china_moral_threat_df,
      family = binomial(link="logit"))
  
  threat_interaction_plot <- rbind(threat_interaction_plot,
                                   data.frame(probs = predict(m_temp, type="response"),
                                              us_greatest_threat = china_moral_threat_df$us_greatest_threat,
                                              neg_morals = china_moral_threat_df$neg_morals,
                                              moralsxthreat = china_moral_threat_df$us_greatest_threat*china_moral_threat_df$neg_morals, 
                                              issue = dvs[i]))
  
}

levels(threat_interaction_plot$issue) <- c("Taiwan Arms Sales", "US Military Presence\nin Asia-Pacific", "US Suppression\nof China", "US Supremacy", "US Support\nfor Tibetan Separatists", "US Recon. near\nChinese Coast")
threat_interaction_plot$us_greatest_threat <- factor(threat_interaction_plot$us_greatest_threat, levels = c(1,0))

# --- Figure A6(B)
ggplot(threat_interaction_plot, aes(x = neg_morals, y = probs, colour = us_greatest_threat)) +
  theme_bw() +
  geom_line(size=1.2) + 
  facet_wrap(~issue, ncol = 3) +
  scale_color_manual(values = c("gray40", "rosybrown2", "lightcyan3"), 
                     name = "US Greatest\nThreat to China?", labels = c("Yes", "No")) +
  labs(x = "Perceptions of American Immorality", y = "Probability of Issue Importance") +
  theme(axis.text.y=element_text(size=12),
        axis.text.x=element_text(size=12),
        axis.title.x = element_text(size=16, margin = margin(t=12)),
        axis.title.y = element_text(size=16, margin = margin(r=12)),
        strip.background = element_rect(fill="gray95"),
        strip.text = element_text(size=14),
        plot.title = element_text(hjust = -.1, size = 26),
        panel.spacing = unit(1.2, "lines"),
        legend.text = element_text(size=14),
        legend.title = element_text(size=14)) 
#ggsave("china_issue_interactions.pdf", width=10, height=7.4)



# --- full tables for issue importance dvs, with immorality x greatest threat as predictor
china_moral_threat_df$education <- factor(china_moral_threat_df$education, levels = c("Less than College", "Vocational College", "Bachelors", "Graduate Education", "No Education Answer"))
threat_interaction_table <- list()
for(i in 1:length(dvs)){
  
  m_temp_simple <-
    glm(eval(as.name(dvs[i])) ~ us_greatest_threat + neg_morals + us_greatest_threat*neg_morals,
        data = china_moral_threat_df,
        family = binomial(link="logit"))

  m_temp_demos <-
    glm(eval(as.name(dvs[i])) ~ us_greatest_threat + neg_morals + us_greatest_threat*neg_morals +
          gender + age + education + income,
        data = china_moral_threat_df,
        family = binomial(link="logit"))
  
  threat_interaction_table[[paste(dvs[i],"_simple",sep="")]] <- m_temp_simple
  threat_interaction_table[[paste(dvs[i],"_demos",sep="")]] <- m_temp_demos

}

# --- Table A2
screenreg(threat_interaction_table)


# --- full tables for issue importance dvs, without the immorality x greatest threat interaction, to show main effects of each
threat_morality_table <- list()
for(i in 1:length(dvs)){
  m_temp_simple <-
    glm(eval(as.name(dvs[i])) ~ us_greatest_threat + neg_morals + 
          gender + age + education + income,
        data = china_moral_threat_df,
        family = binomial(link="logit"))
  threat_morality_table[[paste(dvs[i],"_nointeract",sep="")]] <- m_temp_simple
}

# --- Table A3
screenreg(threat_morality_table)


# --- full tables for issue importance dvs, without the greatest threat predictor, to avoid post-treatment bias issues 
issue_morality_table <- list()
for(i in 1:length(dvs)){
  m_temp_issue <-
    glm(eval(as.name(dvs[i])) ~ neg_morals + 
          gender + age + education + income,
        data = china_moral_threat_df,
        family = binomial(link="logit"))
  issue_morality_table[[paste(dvs[i],"_immoral",sep="")]] <- m_temp_issue
}

# --- Table A4
screenreg(issue_morality_table)



